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Perpetual integral functionals of diffusions and their 

numerical computations 



In this paper we study perpetual integral functionals of diffusions. Our interest 
is focused on cases where such functionals can be expressed as first hitting times for 
some other diffusions. In particular, we generalize the result in j2l] in which one-sided 
functionals of Brownian motion with drift are connected with first hitting times of 
reflecting diffusions. 

Interpretating perpetual integral functionals as hitting times allows us to compute 
numerically their distributions by applying numerical algorithms for hitting times. 
Hereby, we discuss two approaches: 

• numerical inversion of the Laplace transform of the first hitting time, 

• numerical solution of the PDE associated with the distribution function of the 
first hitting time. 

For numerical inversion of Laplace tranforms we have implemented the Euler algo- 
rithm developed by Abate and Whitt. However, perpetuities lead often to diffusions 
for which the explicit forms of the Laplace transforms of first hitting times are not 
available. In such cases, and also otherwise, algorithms for numerical solutions of 
PDE's can be evoked. In particular, we analyze the Kolmogorov PDE of some diffu- 
sions appearing in our work via the Crank-Nicolson scheme. 
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1 Introduction 



Let [Y t : t > 0} be a regular linear diffusion taking values on an interval /. The left and 
right endpoints of the interval are denoted by I and r, respectively. For a locally integrable 
function / : / i— > R + define the perpetual integral functional associated with / and Y via 



f(Y t )dt. (1.1) 

o 

An important example of perpetual integral functionals is 

exp [-2a J dt, a > 0, 

where B^> is a BM with positive drift /x, studied by Dufresne in [TH] in connection with 
risk theory and pension funding. In particular, from [10], this functional is distributed as 
1/(2 a 2 Z v ) where Z u is a gamma-distributed random variable with the density function 

fzA z ) = p^y 2" -1 e"*, i/:=n/a. 

In Yor [HI] (see |32j for an English translation) it is shown 

/ exp {-2aB^)ds = H (R^), (1.2) 
Jo 

where is a Bessel process of dimension 5 = 2(1 — (n/a)) started at l/a, 

H (B®) := inf{t : R? = 0}, 

and — reads "is identical in law with" (in fact, R^ can be constructed in the same 
probability space as B^> and then (ll.2jl holds a.s.). In [21] the methodology used in [SI] 
is developed for more general perpetual functionals for BM with positive drift and, in 
particular, results for one sided functionals are presented. An example of these is 

r 'exp(-2oSW)l {fl M >0} ds = H 1/a (R^), 

J 

where the Bessel diffusion i£CW a ) [ s started at and, in the case < [i < a, reflected at 
0. For further results and references for Dufresne's functionals, see [21], [20], [2H], [21] and 

m 

In this paper, Section 2, we recall (from [I]) the connection between perpetual integral 
functionals and first hitting times. After this, the result in [21], Proposition 2.3, concern- 
ing one-sided perpetual functionals of B^ is generalized for Y (defined via a SDE) and 
functionals of the type in Ijl.lj) . In Section 3, to make the paper more self contained and 
also as an introduction to Section 4, some basic facts about the distributions of the first 
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hitting times are presented. Section 4 contains brief descriptions of the Euler algorithm 
for numerical inversion of Laplace transforms and the Crank-Nicolson scheme for solving 
PDE's, which we implemented in Matlab. The paper is concluded with Section 5 where 
the distributions of some perpetual functionals are computed numerically. In particular, 
we compare the one-sided functionals 



where B t ' := a B t + fit with B the standard Brownian motion. It is also seen that some 
of the diffusions studied have bad singularities making the PDE's numerically troublesome 
to solve. In some cases this problem can, at least partly, be solved by transforming the 
diffusion to a new one with better behaviour. It seems to us that for a general numer- 
ical approach for calculating distributions of perpetualities, more sophisticated PDE or 
other methods such as Monte Carlo simulation are needed for the cases where the Laplace 
transform is not available for numerical inversion. 

2 Perpetual integral functionals as first hitting times 

Consider a diffusion Y on an open interval / = (/,r) determined by the SDE 



where B is a standard Brownian motion defined in a complete probability space J 7 , {J-'t}, P)- 
It is assumed that a and b are continuous and a(x) > for all x e /. The diffusion Y is 
considered up to its explosion (or life) time 



Let / be a positive and continuous function defined on J, and consider for t > the 
integral functional 



Jo 

We remark that {A t : t > 0} is an additive functional of Y in the usual sense (see e.g. [2] 
p. 148). Taking t = ( gives us the perpetual integral functional 




and 




dY t = a(Y t )dB t + b(Y t )dt, 



(2.1) 



C := inf {t : Y t £ I}. 





Assuming that < oo a.s. we are interested in the distribution of A$. 

A sufficient condition for finiteness is clearly that the mean of A^ is finite: 
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where Go denotes the Green kernel of Y and m is the speed measure (for these see, e.g., 
Pj). A neccessary and sufficient condition in the case of a Brownian motion with drift 
fi > is that the function / is integrable at +00 (see Engelbert and Senf [TT] and Salminen 
and Yor |25j). We refer also to a recent paper |22j for such a condition valid for continuous 
/ and a fairly general diffusion Y. 

Next proposition connects the perpetual integral functionals to the first hitting times. 
The result is extracted from Propositions 2.1 and 2.3 in |3] where the proof can be found. 
We remark also that the result generalizes Proposition 2.1 in [2i] . 

Proposition 2.1. Let Y, A, and f be as above and assume that there exists a two times 
continuously differentiate function g such that 

f(x) = (g'(x)a(x)) 2 , xel. (2.2) 

Let {a t : < t < A^} denote the inverse of A, that is, 

at : = min {s : A s > t], te.[Q,A{). 

1. Then the process Z given by 

Z t :=g(Y at ), te[0,A c ), (2.3) 
is a diffusion satisfying the SDE 

dZ t = dB t + G{g-\Z t )) dt, t G [0, A c ). 
where B t is a Brownian motion and 

G(x) = Q o{xf g"{x) + b{x) g'(x^J . (2.4) 

2. Let x & I and y G I be such that P x -a.s. 

HyiY) := inf{t : Y t = y} < 00. 

Then 

A Hy[Y ) = inf{t : Z t = g(y)} =: H g{y) (Z) a.s. 
with Yq — x and Z — g(x). 

3. Suppose g(r) := lim 2 ^ r g(z) exists. Suppose also that the following statements hold 
a.s. 

(i) liml^ = r, (ii) A^ := limAj < 00. 

Then 

A c = H g{r) (Z) a.s. 
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In [2H Proposition 2.3 one sided functional for Brownian motion with positive drift 
are studied. This result is generalized here, under some assumptions, to the present case. 
Suppose G (l,r) and recall that f(x) > for all x G (l,r). Consider the functional 



A° ( := f f{Y s )l {Ys>0} ds. 
Jo 



Let 



C t := / l { Y 3>0 }ds, t < C, 
Jo 

and {c t : < t < B^} denote the inverse of C. We assume also that 

liml^ = r a.s. (2-5) 

It is well known (see [E]) that the process 

Y + := {Y ct :0<t< CJ 

is identical in law with Y living on [0, r) and having as a reflecting boundary point. 
Applying the random time change means that on every sample path the excursions below 
are omitted after which the gaps created are closed by joining the excursions together. 
Therefore, 



A°= [ f(Y+)ds=:A+ 
Jo 



where ( + is the life time of Y + . 

Next introduce the local time of Y + at via 

L t (Y + ) := a 2 (0) lim(2e)- 1 Leb{0 < s < t : Y s + < e}. 

elO 

Under some additional smoothness assumptions on a and b (see McKean [THj) the pair 
(Y + , L{Y + )) with Yq = x > can be viewed as the unique solution of the reflected SDE 

dX t = a{X t ) dB t + b(X t ) dt + dL t (X), X = x, 

such that 

(a) \im t ^ X )X(t) = r, 

(b) < X(t) < r for all t < ((X), 

(c) 1 1 — ► L t (X) is continuous, increasing with L (X) = 0, and 



/ l {0} (X s )dL s (X) = L t (X). 
Jo 
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We now give the promised generalization. 

Proposition 2.2. Let Y + be as given above and define for t < ( + 

At:= f f(Y s + )ds. 
Jo 

The inverse of A + is denoted by {af : < t < A^}. Recall the definition of the function 
g in $2.2}) and define the process Z + via 

Z+:=g(Y+), te[0,A+). (2.6) 

Then 

A+ = inf{t : Z+ = g{r)} a.s. (2.7) 
with Z = g(x). Moreover, Z + satisfies the reflected SDE 

dZ+ = dB t + G(g-\Z+)) dt + dL t (Z + ), t e [0, A+). (2.8) 

where B t is a Brownian motion, 

L t (Z + ) = lim(2e)- 1 Leb{0 < s < t : g(0) < Z+ < g(0) + e}, (2.9) 

ej.0 

and G is as in 42.4)) . The local time L(Z + ) is related to the local time L{Y + ) by 

L t (Z+)=g'(0)L 4 (Y+). (2.10) 
Proof. To fix ideas, assume that g is monotonically increasing. By Ito's formula for u < ( 

g(X+) - g(Y+) = (a(F+) dB s + b(Y+) ds + dL s (Y+)) 

Jo 

+ l£ 9"(Y s + )a 2 (Y+)ds. 

Replacing u by af yields 

Zt - Z+ = £ g'(Y+)a(Y+) dB s + g'(0) L a +(Y + ) 

+ {g'(Y+)o-(Y+)) 2 G(Y s + )ds. 
Jo 

Since af is the inverse of Af and (Af)' = ((7'(F s + )cr(F s + )) 2 we have 

vv = ^ = (^(yjy(yj))" 2 . (2.ii) 
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From Levy's theorem it follows that 

B t := P g'(Y s + )a(Y s + ) dB s , t G [0, A+), 
Jo 

is a (stopped) Brownian motion. Consequently, for t < 

Zt ~ Zt =Bt + g'(0) L 4 (Y+) + f (g'(Y+)a(Y+)) 2 G(Y+) cfa+ 
= B t + £ G{g-\Zt)) ds + g'(0) L at (Y+). 

Clearly, viewing t h- > g'(0) L a +{Y + ) as a functional of Z + then this functional increases 
only on the set {t . : Zf = g(0)}. Moreover, since Y t + > for t > we have Zf > g(0) 
for t > by monotonicity of g. Hence, (Z + ,L(Z + )) can be seen as the unique solution of 
the reflected SDE (TT8l) with L(Z+) as in (fOll satisfying (fOTl . as claimed. Finally, again 
by the monotonicity of g, the identity (|2.7j) follows from the definition (12.61) of Z + and the 
assumption (12 □ 

Remark 2.3. Notice that the above approach yields a stronger result than in i.e., the 
identity \2. 1\) holds a.s. 



3 Reminder on first hitting times 
3.1 Distribution functions and PDEs 

Let Y be a linear diffusion determined via the SDE (|2.1j) . It is here assumed that Y hits 
r a.s. and is killed when this happens. Therefore, the boundary point r is either exit-not- 
entrance or regular with killing, and / is either natural or entrance-not-exit or regular with 
reflection. Letting H r (Y) denote the hitting time of r we have 



/r 
p(t;x,y)m(dy), (3.1) 



where p denotes the symmetric transition density of Y with respect to its speed measure 
m. It is well known (see [IH] p. 149 and [18j) that (t, x) i— > pit; x, y) satisfies for all y G (I, r) 
the PDE 

did 2 d 
—p(t; x, y) = 2 (t2 ( x ) g^P^ x > ^) + x ' ^ 

=: (Gp)(t;x,y) 

and the condition lim x ^ r p(t; x, y) = 0. Moreover, in the case / is regular with reflection or 
entrance-not-exit we impose at / the condition 

d 

hm—p(t;x,y) = 0, 

x^l OX 
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and in the case I is natural the condition 

d 

\imp(t; x, y) = lim —pit; x, y) = 0. 

x^l x^l OX 

See [IH] or [3] for the boundary classification of linear diffusions. 

Letting {T t } denote the semigroup associated with Y we may write from ()3.1|) 

P x (H r (Y)>t) = T t l(x). 

Recall from [TBI (where the case with natural scale is treated) that (t,x) i— > T t g(x) with g 
bounded and continuous satisfies 

^ j {T i g){x) = {QT t )g{x). 
Consequently, the distribution function 

(t,x)»u(t,x) :=P X (H Z (Y) <t) 
is the unique solution of the PDE problem 

^u(t,x) = (Qu)(t,x) (3.2) 

with the initial condition \im t ^ou(t,x) = for all x G (l,r) and the boundary condition 
lim^,, u(t, x) — 1 for all t > 0. Further, if I is regular with reflection or entrance-not-exit 

d , v 
lim —u[t, x) = 0, 

x^l OX 



and in the case I is natural 
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\vmu(t,x) = lim— — tt(t, x) = 0. 

x->/ x-t-J OX 

Remark 3.1. Using the fact (see J75j/j i/iai 

i-> ^p(t;x,y) 

is continuous and satisfies the same boundary conditions as the density p(t; x, y) it is easy 
(at least when m(l,r) < oo) to deduce that 

|p.(jr.(y)<«) = -jh 55 ^A p( t;,,v), 

where 

S'(y) = exp (- J y 2 o-- 2 (v) b(v) dv 
is the derivative of the scale function S (cf. fTb^ p. 154). 
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3.2 Laplace transforms and ODEs 

For the approach with the Laplace transform of H r (Y) consider the second order ODE 

Qu{x) = Xu(x), (3.3) 

where A > 0. It is known (see [12] p. 488, and [HU p. 128) that the equation (13. 3 j) has a 
positive increasing solution ip\ and a positive decreasing solution (p\. In case I is natural 
or entrance and r is exit these solutions are unique up to multiplicative constants. When 
I is regular with reflection the condition V'aCO = must be posed, and when r is regular 
with killing the condition is <px(r—) = 0. The Green kernel G\ of Y can be expressed via 
these solutions as 

poo 

G x {x,y) := / e~~ xt p{t;x,y)dt 



'0 



-f- i/)x(x)tpx(y), x<y, 



where w\ is the Wronskian (see e.g. Using the Green kernel the Laplace transform 

for the first hitting time H y {Y) is given by 



and, in particular, 



E, (e-™M) = ^) . (3.4) 



4 Numerical methods 

4.1 Numerical inversion of Laplace transforms 

There are several efficient methods for numerical inversion of the Laplace transforms of 
probability density functions or probability distribution functions. We have implemented 
a method developed by Abate and Whitt in [T]. This so called the Euler-algorithm has 
proved to be very effective in many applications see e.g. [12], pj, and |Hj. The main 
features of this method are presented below. For more details and also for further references, 
see [l]. 

Consider a non-negative random variable with density / and its Laplace transform 

/■oo 

/(A):= / e- xt f(t)dt. 
Jo 

The well known inversion integral formula (called the Bromwich or also the Fourier-Mellin 
integral) states that 

-r />a+ioo 

f(t) = — e xt f(X)d\, (4.1) 



2?ri 



a— l oo 
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where it is assumed that / does not have singularities on or to the right of the vertical line 
A = a. 



Remark 4.1. For first hitting times of diffusions considered in Section 3 we have (cf. 

/(A) = E x (e-**W) = ^4. 

If the left boundary point I is not natural it follows from the classical theory of second order 
differential operators (see e.g. thatip\(x) is for every x G (l,r) an entire function of X 

and the zeroes o/A h» ip\(x) are for every x G (I, r) simple and negative. Consequently, the 
inversion formula U-l]) holds in this case for any a > 0. If I is natural we can approximate 
the first hitting time H r (Y) via a sequence of first hitting times {H r (Y^ n ')} associated with 
the diffusions Y^ n \ n = 1,2, ... , constructed from Y by reflection at I + K respectively. 
Then H r (Y ( - n ^) — > H r (Y) in distribution and by dominated convergence it is seen that \4.1\) 
is valid also in this case. 

Since Re(/(a + iu)) = Re(/(a — iu)), Im(/(a + iu)) = — Im(/(a — iu)), and f(t) = 
for t < the inversion integral (|4.ip takes the form 



f{t) = / Re(/(a + iu)) cos{ut)du. (4.2) 

7T Jo 



Next we approximate f(t) by using the trapezoidal rule for the integral (|4.2jl (see [T] for 
some comments on the effectiveness of this procedure). Letting h denote the step size we 
have for fixed a and t 

hp at - ?he at J 30 , 

m*h(t) = Re(/(a)) + — Re(/(a + i)) cos(Ef). 

k=l 

Choosing h = n/(2t), a = A/(2t) (with A to be made precise later) and truncating the 
infinite series to the first j terms we are led to define 

p a / 2 

EC" 1 )* *®' ( 43 ) 

fc=0 

where 

a (t) :=f(A/2t) /2 

and 

a k {t) :=Re(/((A + 2fc7ri)/2t)) , fc = 1, 2, . . . , j. 

It is possible to accelerate the convergence of the series in ((Qj) by considering it as an 
alternating series (which is not the case in general) and using the Euler summation with 



10 



binomial weights (see p]). Hence, the proposed final approximation with parameters m, 
n, and A is 

fit) « £(m,n,t) := X) (T) 2_W W*), 
fc=0 ^ ' 

In the examples below we use, following pQ, m = 11 and n = 15. The error associated with 
Euler summation can be estimated by considering the difference 

E(m, n + 1, t) — E(m, n, t). 

It is advantageous from numerical computational point of view to invert, instead of the 
density function, the complementary distribution function. Therefore consider 

™<x> 



POO 

F C (A): = / e- xt (l-F(t))dt, 
Jo 



where F is the distribution function associated with /. Firstly, the fact |1 — F(t)\ < 1 can 
be used to show (see fl)) that 



where stands for the discretization error when approximating the integral in ()4.2|) for 
F c via the trapezoidal rule. For instance, A = 18.4 gives the upper bound 10" 8 . Secondly, 
under some additional smoothness assumption, it can be proved (see [I] Remark 1) that for 
F c we have a^it) > when k/t is large enough motivating the use of the Euler summation 
(since the series in (|4.3jl is now alternating). 

For the first hitting time of r for the diffusion Y the Laplace transform of the comple- 
mentary distribution function is given by 

nX) = 1 (1 - E, (e— 00)) = 1 _ 
_ 1 ipxjr) - j)\{x) 



4.2 Numerical solutions of PDEs 

In this section we describe using [2H] and [Uj two finite difference methods known as the 
Crank- Nicols on (C-N) scheme and the backward Euler (BE) method. The C-N scheme is 
used with satisfactory results in [21] for calculation of transition probability densities of 
certain diffusions. The BE method can be applied in connection with the C-N scheme for 
the first time step to damp some numerical oscillations typical to the C-N scheme. Both 
methods are unconditionally stable: the numerical solutions are well behaved (do not blow 
up) for any choice of At. Because of the singularities appearing in the drift coefficients of 
our equations, methods that do not have this property (such as the explicit, forward Euler 
method) are practically unusable since they would require extremely small time steps. 
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To start with, let us introduce a uniformly spaced grid on the rectangle [0, T] x [l,r] 
with (M + 1) x (TV + 1) nodes, that is, for 

T \ r ~ l 
At = — , Ax = , 

we let t m = mAt and x n = 1 + nAx where m — 0, 1, M, n — 0, 1, N. For a real valued 
function / we use the following finite difference approximations, which can be justified with 
Taylor's expansion: 

• the forward difference approximation of the derivative of / is 

■ the centralized difference approximation of the derivative of / is 

| (l() = Z(£«lzi(£tl) +0((Al)2) . 

• the centralized difference approximation of the second order derivative of / is 

The C-N scheme approximates the left hand side in equation (|3.2jl with the forward dif- 
ference and the right hand side with the average of the centralized differences at two 
consequtive times. Denoting u™ := u(t m ,x n ) and dropping the truncation error terms, the 
discretized equation then reads 

M n' +1 ~ U ri \ 2( \ \ ( U n+l ~ + U n ^ ~ + U™_ 1 

At 2 a{Xn) 2\ (Ax) 2 + (Ax) 2 

1 /V m + 1 — V m + 1 a,™ — V m 



+ b(x n )-( n+1 n ~ l + 

v ! 2\ 2Ax 2Ax 

Multiply both sides with At, and define r\ = r 2 = 2 (aI) 2 m R earran g m g the terms so 
that the values at time t m+ \ appear on the left hand side and values at time t m appear on 
the right hand side, we have 

where 

^ = ^(b{x n )ri -a 2 (x n )r 2 ), (4.4) 

5 n = l + a 2 (x n )r 2 , (4.5) 

Cn = ^(b(x n )n + <x 2 (x n )r 2 ), (4.6) 

D n = l-a 2 (x n )r 2 . (4.7) 
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Together with the boundary conditions, these form a set of N + 1 linear equations which 
we then solve for each m = 1,2, ...,M + 1, using the initial condition for m — 0. The 
Neumann boundary condition (which is needed at a reflecting or an entrance boundary 
point) is implemented with the second order approximation. In other words, from 

= 



2Ax 

we have = U™ for the value of u at a "ghost" point beyond the boundary. Plugging 
this into the discretized equation at n = gives 

(1 + r 2 a 2 (x n ))u™ +1 - r 2 a 2 (x n )u? +1 = (1 - r 2 a 2 (x n ))u™ + r 2 a 2 (x n )u?. 

Using a second order approximation for the boundary condition seems to be important in 
order to have good convergence. 

For the backward Euler method, one similarly takes the central approximations for the 
spatial derivatives but now only at the time step m + 1. This leads to the equation 

A 7,™+! _|_ R _ n 7 ,m+l _ „m 

for n = 1,...,N and m = 1, ...,M, where A n , B ni C n are given in (|4.4jl - ()4.6jl but now 
with T\ = r 2 = , A \ 2 . The second order implementation of the Neumann boundary 
condition becomes 

(1 + r 2 a 2 (x n )X +1 - r 2 a 2 (x n )u? +1 = v%. 

Both methods described above are second order accurate in space, but only the C-N 
scheme is second order accurate in time. However, C-N is known to produce numerical 
oscillations around discontinuities and sharp gradients if the drift term is large compared to 
the diffusion coefficient (see for example [H], [Hj), while the BE method does not have this 
problem. In our examples this is seen as oscillations near the killing boundary. Although 
these oscillations were damped quite rapidly both in space and time, they still produced a 
small phase shift in the numerical solution of the hitting time distribution t i— ► u(x,t). For 
this reason the first step in the C-N scheme is divided into 10 — 100 substeps, as suggested 
in [H]. For the substeps we used the BE method. While small oscillations still remained 
for some of the examples, this procedure provided sufficient damping to give very accurate 
results in our test cases. 

A problem more serious than the oscillations appearing in the C-N scheme is faced 
when the diffusing particle is pushed away from an exit boundary by a drift term tending 
to infinity in the vicinity of the boundary. In such cases convergence can be very slow 
or even nonattainable without huge computer capacity. We discuss this problem in more 
detail in the examples below. 



5 Examples 

In this section some perpetuities are examined numerically. If the Laplace transforms of 
the functionals are available we use the Euler algorithm for computing the density and/or 
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the distribution functions. Applying in these cases also the numerical methods based on 
the associated PDE's obtained from the hitting time representations of the functionals and 
checking that the solutions resulting from the different methods coincide we are able to 
verify the correctness of the implementations. 

Numerical methods for solving PDEs constitute a powerful tool for computing hitting 
time distributions of diffusions in general. However, there are diffusions as in Example 15.41 
for which the methods presented here work unsatisfactorily. At least in some particular 
cases it is possible to transform the diffusion to a new one for which the methods seem to 
work better. This is discussed in Example 15.41 Due to these difficulties one may consider 
the numerical inversion of the Laplace transforms when available as the first choice for the 
kind of numerical computations studied in the paper. 

Below the notation {B[^ : t > 0} is used for a Brownian motion with (infinitesimal) 
drift fi and variance <r, i.e., = a B t + fit with B a standard Brownian motion. We 

assume that a > and fi > 0. In case a = 1 we write for B^ 1 ' 1 '. For a Bessel process 
with dimension parameter 8 we use the notation of {Itf : t > 0}, and refer to [H] for their 
properties. If nothing else is written it is assumed that B^ 1 '^ = and = 0. 

Example 5.1. Consider the perpetual integral functional 



The Laplace transform of this functional is computed in |3] and f° r an arbitrary initial 
value x; taking therein x = gives 





where 



and 2-F\ denotes Gauss's 




OG 



r(a + k) r(6 + k) x k 

r( c + k) T\ 



2 F 1 (a, b, c; x) : 



T(a)T(b) 



E 



a(a + 1) . . . (a + k - 1) b(b + 1) . . . {b + k - 1) x k 
c(c+l)...(c + fc- 1) ~k\' 



k=l 



Applying Proposition 12.11 with g(x) := 2arctane a: we obtain 



h = H n (Z) 



a.s., 
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where Z satisfies 



dZ t = dB t +[\ ctn Z t + — ^— ) dt, Z = tt/2. (5.1) 



2 sin Z t 

For the drift term in l|5.1jh it holds 

^ -\t ss 1 V 1 ( l^ a; 1 / 1\ Z 

Gig (x)) = - ctnx H = - u tan — I — u H — ctn — . 

vy v 2 sinx 2 V 2j 2 2V 2/ 2 

Notice that for < /i < 1/2 the drift of Z tends to — oo when Z is approaching n. 

In Figures l5~Tl and HI we present the density and the distribution functions, respectively, 
of I\ computed with the Euler algorithm. It has been checked that the PDE method yields 
the same results. However, the convergence in the case \i = 0.3 seems to be slow. In fact, 
for n — 0.1 the convergence rate of the PDE method is so slow that we were unable to get 
satisfactory results with our limited computation capacity (RAM). 



Example 5.2. In this example we compare the functionals 

/"OO /*00 

h ■= / exp(-2B^ a) )ds and J 3 := / (exp(B^ a) ) + l)~ 2 ds. 
Jo Jo 



Notice that 



J 3 = / exp(-2fi^' CT) ) — ds, 

Jo (1 + exp(— 



exp(- J B^' CT) )) 2 

hence I3 may be seen as a modification of I2 which does not allow arbitrary large positive 
discounting. We remark also that I 3 has all moments which is not the case with I 2 - 
The Dufresne-Yor identity (cf. Q1.2|) ) states that 

h = H (R> 2 - 2 ^) a.s., 

where R 2 2 ^ CT ^ = 1/cr. Consequently, 

M«p(-p*)) = 

with (see |2j p. 133) 

<p p (x) = x~ v K v (xy/2p), and <^ p (0) = 2^ +2 )/ 2 r(-z/) p^ 2 

and z/ = —fi/a 2 . 

For 7 3 we have the identity 

J 3 = flo(Z) a.s., 
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0.45 




5 10 15 



Figure 1: Density of J °° cosh 2 {Bf'')dt for \x = 0.3 (lowest peak), /x = 0.5 and 
(highest peak). 




Figure 2: Distribution function of J °° cosh 2 (B t ) dt for \x = 0.3 (lowest curve), 
and /j, = 0.7 (highest curve). 
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where Z is the diffusion associated with the SDE 

1 exp(Z t 



dZ t — dB t + I p + (p — -)- ^-^t) dt, Z = -\og2. 

Notice that here g(x) := — log(l + exp(— x)) (cf. Proposition EH]) and that Z lives on R„. 
From |lj we recall the Laplace transform 

E (exp( - p/ 3 )) = K 2 F 1 (a,(3, a + (3 + 2p ; 1/2), 

where 



a=^-p+ v / /i 2 + 2p+yi + 2p, /3 = I-p + vV + 2p-yi + 2p, 

and 

_ r(2p + q) T(2p + /j) 
T(2p + a + (3) T(2p) ' 

See Figures EJ IU E] and El for illustrations of the distributions of I 2 and I 3 computed with 
the Euler algorithm. 

For both functionals in this example it was possible to solve the corresponding PDE 
numerically for p > |. For p < 1/2 the drift term tends to — oo as Z approaches the 
killing boundary 0. This again leads to very slow convergence. While it was still possible 
to achieve good results for some choices of p < ~, for a small enough p the results were 
bad even with the finest grid we could run on the computer. Notice that here we also 
need to truncate the semi-infinite domains into finite ones for numerical computations. 
This did not constitute a major problem, but with larger domains it is difficult to achieve 
(depending on the computer capacity) a grid which is spatially dense enough for accurate 
computations. 

Example 5.3. We define the one-sided variants of 7 2 and I 3 via 



POO 

I A := / exp(-25^ <T) )l 
Jo 



and 



POO 

J (exp(i^)) + l)-H {B ^ )>0} ds, 



respectively. 

In [21] it is shown that 

h = H 1/a (R^ /a2) ) a.s. 
where ^ = 0. The Laplace transform of h is hence given by 



E (exp(-pL 
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0.5 1 1.5 2 2.5 3 

Figure 3: Density of J" °° exp(— 2 b\^) dt for (i = 0.75 (lowest peak), fi = 1.50 and fi = 2.50 
(highest peak). 




Figure 4: Distribution function of J °° exp(— 2 Bf"') dt for /x = 0.75 (lowest curve), /i = 1.50 
and ji = 2.50 (highest curve). 
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2.5 




0.5 1 1.5 2 2.5 3 3.5 4 



Figure 5: Density function of J °°(exp( B\ f') + 1) 2 dt for fi = 0.25 (lowest peak), // = 0.50 
and /i = 0.75 (highest peak). 




Figure 6: Distribution function of J °°(exp(5 t ) + 1) 2 dt for fi = 0.25 (lowest curve), 
H = 0.50 and = 0.75 (highest curve). 
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0.25 



0.15 




Figure 7: Density function of f£°(exp(B. 
with the density function of / °°exp(— 2B 



(/i,cr)x 



\-2 



1 {B< t ^ a '>>0} 



1 {B t ( "' CT) >0} 



dt (upper peak) compared 
dt for n = 0.04 and a = 0.20. 



with (see [2J p. 133) 

ip p (x) = x~ v I v (xy/2p) and V/>(°) 



P 



y/2 



and v = fi/cr 2 — 1. 

The Laplace transform of the functional Is (in [21] this is called the one-sided translated 
Dufresne functional) is not known but the following identity (see [23]) holds 

J 5 = H Q {Z) a.s., 

where Z is a diffusion associated with the generator 



yf(x) = ~ — (x) + ( -a 



2dx 2 



2 a (1 — exp (ax)) J dx 



df, 



x 



living on [—(log 2)/cr, 0), having — (log2)/cx as a reflecting barrier, and as a killing barrier. 

In Figure Owe compare the densities of 7 4 and I 5 (see |9j for comparisions between I 2 
and 7 4 ). The density and distribution functions of 7 5 are displayed for different values of 
\i and a in Figures El and El 

Example 5.4. In our final example we consider the functional 

POO 

I ( 6 5) := / exp(-2 RW)ds, 5 > 2. 
Jo 



20 



0.35 




Figure 9: Distribution function of J °°(exp(5 t ) + 1) 2 lr B M >0 j dt for a = 0.20 and 
\i = 0.03 (lowest curve), fi = 0.04, and /x = 0.05 (highest curve) 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 



Figure 10: Density function of J °° exp(— 2 R\ ') dt for 5 = 3.0 (upper peak) obtained using 
the drift ^(1 + compared with the correct one (lower peak). 

Proposition 12. II when applied for and g(x) := exp(rr) leads us to the identity 

h = H (Z) a.s. (5.2) 
with Z a diffusion associated with the SDE 

dZ t = dB t + -^-(l + dt, Z = 1. (5.3) 

2z t v log^y 

In the case 5 = 3 it is known (see Legall [TJj and also [21]) that 

4 3) = H^R®) a.s. (5.4) 

with R® = 0. 

Since we do not have an expression for the Laplace transform of 1^ for 5 ^ 3 we solve 
numerically the associated PDE. Unfortunately, due to the complexity of the drift term 
(in particular, notice that this tends, for all values on 6 > 2, to +oo in the vicinity of 0+) 
simple finite difference schemes do not seem to give solutions converging to the correct one, 
see Figure [ni In search for improvement we implemented a nonuniform grid making the 
spatial discretization denser near the boundaries, and used a fourth-order implementation 
at the Neumann boundary. While this yielded better results than what is seen in Figure 
ITTH full convergence still remained out of reach. 

These difficulties can at least partly be overcome by transforming the diffusion Z given 
via SDE (|5.3jl . Indeed, we study now the /i-transform of Z with h(x) = S(x) — S(0) where 
S is the scale function of Z Straightforward computations (cf. [3j p. 17) show that we may 
take | 

S(x) = |logx| 2_<5 , < x < 1 

o 2 
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(for simplicity we consider only the case 5 > 2) . Then 



lim S(x) = and S'(x) = x 1 | logo: 



is 



Consequently, the generator of the /i-transform is given by 



id 2 f 




5-1 
logx 
3-5 
logx 



df , S'(x) df 



2 dx 2 
Id 2 / 



dx S(x) dx 



2 dx 2 



df 

-f, 0<x<l. 



Let denote the /i-transform, i.e., is the diffusion associated with the generator . 
By Williams [30j time reversal result (see |Sj p. 35, also for further references) 



The PDE associated with Z^ seems to be well suited for numerical computations. 
Notice, in particular, that if 5 > 3 the drift term of Z^ tends to +oo as x — > 1 which fact 
is in strong contrast with the corresponding behaviour of the drift term of Z. Hereby it is 
also of interest to classify the boundaries of Z and Z\ It holds for Z that the boundary 
point is exit-not-entrance and 1 is entrance-not-exit. For the process Z^ we have that 
is entrance-not-exit and 1 is entrance-exit (regular) if 2 < 5 < 4 and entrance-not-exit if 
5 > 4. Figures HU H21 show the density and distribution functions of Iq for some choises of 
5 computed from the PDE associated with Z^ . 

As a final comment, and as an extra bonus from our transformation, we remark that 
when 5 = 3 then is the generator of R^ 2 \ and we have recovered the identity (|5.4|) . 

Acknowledgements. Paavo Salminen thanks Vadim Linetsky for information on 
numerical inversion of Laplace transforms of hitting times. Olli Wallin thanks Siddhartha 
Mishra for several useful discussions on numerical solution of partial differential equations. 
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0.5 1 1.5 2 2.5 3 




Figure 12: Distribution function of J °° exp(— 2 R t ) dt for 5 = 2.5 (lowest curve), 5 = 3.0, 
5 = 3.5 (highest curve). 
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